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1 Introduction 


Inverse scattering problems, originally introduced for potential scattering in quantum physics (see for instance Ref. jt] 
and the references therein), have found their applications in a wide variety of scientific and technological fields. An 
important limitation to the application of inverse scattering methods to real data is related to the inherently unstable 
behavior of the algorithms involved [^, [|. We present in this paper a novel method to solve the one-dimensional 
inverse scattering problem, and moreover a discussion about the stability of the obtained solutions is given. 

The inverse problem for the one-dimensional Schrodinger equation can be solved by using the Marchenko equation. 
Burridge showed that the inverse problem of the one-dimensional wave equation can be transformed into an quantum 
inverse scattering representation jj]. This result made it possible to solve electro-magnetic, optical and acoustic inverse 
problems by using the Marchenko equation. 

The Marchenko equation maps a reflection time series, (hereafter to be referred to as the “data-set”), onto an 
integral kernel from which the quantum mechanical potential function can be recovered. The data-set is the Fourier 
transform of the reflection coefficient. Sabatier has discovered a method to solve this integral kernel in the Fourier 
domain for reflection coefficients that are a rational function of the wave number J5|. 

In Ref. H im stability estimates are made for the solutions derived by Sabatier. A disadvantage of the method 
presented in Ref. Q is that the method is highly dimensional and therefore not easily implementable for realistic 
data-sets. In this paper we re-open the discussion with respect to the solvability and stability of one-dimensional 
Marchenko-equation on the basis of novel results which have recently become available in the field of the integrability 
of nonlinear evolution equations || 10 . It is shown in these publications that large classes of solutions of integrable 
nonlinear evolution equations can be obtained by solving the Fourier expansion coefficients of its solution by using 
a recursion relationship. Since it is well known that large classes of nonlinear evolution equations are integrable by 
inverse scattering transformations, it is worth investigating whether the machinery developed in Ref. [0, [l(J can also 
be applied in solving quantum inverse scattering problems. We show in this paper that it is possible to transform 
the Gelfand-Levitan-Marchenko equation into a linear recursion relationship which determines the Fourier expansion 
coefficients of the kernel. This result made it also possible to investigate the stability of the solutions by using Lyapunov 
exponents which can be computed from the data-uncertainties. This implies that we have a simple tool which requires 
only information about the data and data-uncertainties, to investigate whether the Marchenko equation is stable. 
Secondly, if we might conclude that the inverse scattering problem is unstable, we provide machinery to extract the 
largest sub-set from the scattering data which leads to a stable solution. 

This paper has the following structure: In Section 2 we focus on the linearization method which leads to the 
relationship between the expansion coefficients of the data and the expansion coefficients of the kernel. We firstly 
investigate the principle for a data-set having one Fourier component only, but the method is generalized to data¬ 
sets having a arbitrary (but finite) number of Fourier components. In Section 3, we investigate the stability of the 
obtained solutions by using Lyapunov exponents. A simple numerical example is given. The paper is concluded with 
a discussion. 


2 The one-dimensional quantum inverse scattering problem 

In this paper, solutions of the quantum inverse scattering problem are discussed. The quantum inverse scattering 
problem is the inverse problem of the Schrodinger equation which is given by: 

tp"(k , x) + k 2 ip(k , x) = V (x)'i p(k, x) (1) 

For reasons of simplicity the following restrictions on the scattering solutions are imposed. It is assumed that V ( x ) : 
1R + —> M , and secondly, only incoming waves from the left are taken into account. The solutions of the Schrodinger 
equation ( 0 ) satisfy: 


ip(k,x) 


e ikx + R(k)e~ ikx x —» —oo 
T(k)e' lkx x —> +oo 


( 2 ) 
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where T(k) is the transmission coefficient, and R(k) is the reflection coefficient from the left. In the following, the 
data-set is defined by the function A(t). The relation between the data-set A(t) and the spectral reflection coefficient 
R(k) of the Schrodinger equation is given by: 


/ OO 

dk R(k)e lkt + bound states 

-OO 

The inverse problem for the Schrodinger equation can be solved using the Marchenko equation: 

/»00 

K(x,y) + A(x + y) + / dz K{x, z)A(y + z) = 0 

J X 

The potential is recovered from the integral kernel K(x,y) by [jllf : 

V(x) = —2—K(x,x) 
dx 


( 3 ) 


( 4 ) 


( 5 ) 


In this paper, we firstly discus a method to obtain solutions of the Marchenko-equation. The method is based upon 
results developed in Ref. |g| . As a starting point, we assume that the data A(x) can be expanded in a Fourier series: 


N 


A{x) = y, .4,1 


ik n x 


( 6 ) 


By considering a data-set which can be expressed in the form as presented in Eq.(jGj), we have implicitly assumed that 
the reflection coefficient is a rational function of the wave-number. This is not a limitation for the application to real 
data, since a measured data-set can be approximated sufficiently by considering a series of the form (^|) for which 
N —> oo. As a first step, we consider the most simple case in which the data-set A{x) has only one Fourier component: 


A[x) = Ae ikx 

In this special case, we propose a kernel K[x,y) having the following structure: 

OO 

K(x, y) = Y B(n,m)e inkx+imky 

n,m =1 

If the expressions (0) and (||) are substituted in Eq.([|), we obtain the following result: 

oo oo /»oo 

Y B(n,rn)e ik( - nx+my) + Ae ik{ - x+v) = - Y AB ( n ’ m)e inkx e ikv / dze i{m+1)kz 

„ -_1 J X 


( 7 ) 


( 8 ) 


( 9 ) 


n,m= 1 


n,m=l 


If we assume that A: is a complex number with positive imaginary part, the integral in Eq.(^|) can be replaced by the 
primitive of an exponential function: 


E B(n m \ e ik ( nx+rnv '> + Ae ik( ' x+V ' 1 = AB(n , m) ^j( n+m+ i)kx c iky 

y i( m +l)k 

n,m=l 


( 10 ) 


71 , 171=1 


We can solve Eq. (0) by comparing the exponential functions of equal powers on both sides. This is done by firstly 
defining a number M = n + m, and hence by comparing all the exponential functions for which M = 2,3,4, ••• 
respectively. If we use the fact that the coefficients A and k are determined by the data-spectrum, the kernel K(x,y) 
can be computed by collecting all the coefficients B(n, to). From the y-dependence on the right-hand side of Eq. 0), 
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we can immediately conclude that all the coefficients B(n,m) vanish for which m ^ 1. We consider firstly the case 
that M = 2: 


Ae ik{x+y ) + B{ 1, l)e ife(x+s/) = 0 


( 11 ) 


This leads to the result that B(l, 1) = —A. If we put M = 3, it can easily be verified by substitution that both the 
coefficients B(l, 2) and B( 2,1) vanish. In general, only coefficients B(n , m) for which m = 1 are unequal to zero and 
follow from the following recursion relationship: 


B{n + 2,1) = A 


B{n, 1) 
2ik 


By using Eq.([f|), we can compute the kernel B(n , 1) as a function of A and k: 


K(x, y) = - Ae ik{x+v) + 


D 5 ikx-\-iky 


2 ik {2ik) 2 


( 12 ) 


(13) 


The potential function V(:r) can be recovered from the kernel K(x,y) by using the relationship (||). If we substitute 
k = i/3 and A = 2d into Eq.(|l3|), we obtain the following expression for K(x, x ): 


o rl 2 9W 3 

K(x, x) = 2de~ 20x + r|_ e -4/3x + - e -6/3x + . . . ( 14 ) 

t> p 

The solution of the kernel (0) is in agreement with solutions obtained by using Sabatiers approach II- The 
advantage of the method presented in this paper is that we solve the inverse scattering problem by deriving a linear 
recursion relationship between the Fourier expansion coefficients of the data and the expansion coefficients of the 
kernel. This implies that we have discovered a method which is easy to implement from a computational point of 
view, and moreover, the stability of the obtained solutions can be investigated by using Lyapunov exponents. 

In order to generalize this method to the data-set ( 0 ), it is also illustrating to discus a data-set consisting of two 
Fourier components: 


A(x) = A ie lklX + A 2 e ik2X 


(15) 


In this special case we propose a kernel K(x,y) having the following structure: 


K{x,y)= J2 B(p,q,r,s)e ipklx e iqkiy e irk2X e isk2y 

p,q,r,s= 1 

If we insert Eq.(|l5|) and Eq.(|l6|) in the Marchenko equation, we obtain the following result: 


OO 

Aie ikAx+y) + A2e ikA*+y) + B( K p,q 1 r,s)e ipklX e iqklV e irk2X 

p,q,r } s —0 


e isk 2 y 


Ai B(p,q,r, s ) 
i([q+ l]ki + sk 2 ) 

m p,q,r,s =0 VL J 7 


giip+q+^kxx e i(r-+s)k 2 x gikiy 


A 2 B(p,q,r,s) 
^ n*(I s + l ] fc 2 + qki) 

p,q,r,s —0 


e Hp+q)k 1 x e i(r+s+l)k 2 x e ik 2 y 


(16) 


(17) 


If we use the fact that the coefficients Ai and A 2 are determined by the data, we can compute the non-zero coefficients 
B(p , q, r, s). This is done, by firstly defining M = p + q + r + s and comparing the exponential functions for which 
M = 2,3, • • • respectively. If we consider the case that M = 2, we obtain the following result: 

Aie ikAx+y) + Aae ik 2 (x+y) + Q^ik^x+y) + ^ ^ X yik 2 (x+y) = Q ( 18 ) 
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From Eg. (jl§|) . we obtain that 5(1, 1,0,0) = —A\ and that 5(0,0,1,1) = — A 2 . In general, it follows from Eq. (0) 
that all the expansion coefficients 5(p, q, r, s) can be computed by solving the following linear iteration relationship: 

B(p,q,r , s) 


B(p + q + 1,1, r + s, 0) = A 1 


B(p + q, 0,7- + s + 1,1) = A 2 - 


i([q + l]fci + sk 2 ) 
B(p , q,r, s) 


(19) 


'*([« + 1 ]k 2 + qk\) 

It can be witnessed from Eq. © that the number of expansion coefficients double after each iteration. As an example, 
we compute the coefficients arising from the first three iterations: 


5(3,1, 0,0) = Ai B(1|1 ’ 0 ’ 0) 


2iiki 


5(1,1,0,0) = -Ai 
5(0,0,1,1) = -^ 2 


5(l,l,2,0) = A 1 a( ° 2 f fc ’ i 1 ’ 1) 

5(2,0,1,1) = 1 2 M 

5(0,0,3,1) = A 2 B(0 2 ^ 2 1 - 1) 


5 ( 5 , 1 , 0 , 0 )= A, B ^ fc ’ i °' 0) 
5 ( 3 ,l, 2 , 0 )=A 1 B( 1 2 t 1 fc ’ i 2 ’° ) 
B (3,1-2,0 

5(1,1,4,0) = ^gg 
5(4,0,l,l)=A 2 Mf 

5(2,0,3,l) = A 2 ^|f 

5(2,0,3,l) = A 2 ^|^i) 

5(0,0,5,1) = 1 2 M 


( 20 ) 


It can be witnessed from Eq. (|20|), that some of the coefficients B(p , g, r, s) occur more than once. Since these coefficients 
all belong to the same exponential basis function, we can simply add them together. Similarly as for the case that the 
data-spectrum has only one Fourier component, we have a iteration relationship to solve the coefficients 5(p, q, r, s). 
If we collect all the coefficients B{p, q , r, s), we find that the kernel K(x, y) is given by the following relationship: 

2 2 . . 2 


K(x,y) = A i^ 


ki(x+y) _ 


*£ 

i,j =1 


A,A 


h 


i_j_ ikj(x+y) 2ikiX _ 


1 E 

i,j,l =1 


AiAjAi 


(kj + ki)(ki + kj ) 


e iki(x+y) e 2 i(k i +k j )x _|_ . . . ( 21 ) 


This result can be made general by considering the full data-set A(x) as presented in Eq.(j(jj). In this case we propose 
a kernel K(x, y) having the following structure: 


Q in\k\x ^h\k\z 


_ giriN kw x ginjv -2 


If we substitute Eq. 

N 


K(x, y) = ^ 5(ni,hi • • •nw,h-Ar)e l: 

ni,ni”-niVj^JV=l 

and Eq.([|) into the Marchenko equation we obtain the following expression: 


( 22 ) 


J2 A n e ik " {x+v) + J2 B(ni,fii • • ■ n N , h N )e iniklX e ifllklV ■ 


, ginNkNXgifiNkNU 


,njv = l 
00 


£ 


A 1 B(n 1 ,h 1 ---n N ,h N ) ^ 

+ l]fei + h 2 k 2 + • • • hjvfc/v) 


i(m-\-hi-\-l)kix i(n2+fi2)k2X 


. e i(riN+n N )kNXgikxy 


A N B{m,h\ ■ • • n N ,n N ) 

i(n\k\ + n 2 fc 2 + • • • [tin + l]fcjv) 


e i(n 1 +n 1 )k 1 x e i(n2+n 2 )k2X i 


, e i(n N +n N + l)k N x £ ik N y 


From this result, we can immediately derive that: 

5(1,1,0,0, • • • , 0,0) = -A 1 
B{ 0,0,1,1, • , 0,0) = — A 2 


(23) 


(24) 


5(0,0,0,0, • • • , 1,1) = —A 


N 
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Moreover, we can derive by using a similar argumentation as used in the case that N = 2 that the recursion relations 
are given by: 


B(ni +h i + l, 1, n 2 + n 2 , 0, ■ • • ,n N + h N , 0) = A 
B(ni + 77,1 ,0, n 2 + h 2 + 1,1, • • • ,n N + ii N , 0) = A 


n 2 , n 2 , ■ ■ ■ ,n N ,n N ) 

1 + l]jfei + fi 2 k 2 + ■ ■ ■fiN^N ) 

B(ni,hi,n 2 , h 2 , ■ • • , n N , h N ) 

i{h\k\ + [n 2 + 1 }k 2 + ■ ■ -n N kN) 


(25) 


B(n i + hi, 0, n 2 + h 2 , 0, • ■ • 


, tin + tin + 1 , 1 ) = A n 


B(ni,hi, n 2 , n 2 , ■ • • , n N , h N ) 
i(fiiki + h 2 k 2 + ■ • • [iijv + l]fc/v) 


Once the coefficients B{n\, n±, n 2 , n 2 , ■ ■ ■ ,riN,nN ) are determined, the kernel K(x,y) can be uniquely recovered. 
Hence the kernel K(x,y) can be represented in the following form: 


N 

K{x,y) =Y J A i e iki ^ +y '> + 

i= 1 


N 

^k, 

*,7 = 1 


A,A 


3_ e ikj(x+y) e 2ikiX _|_ 


N 

E 

1 


AiAjAi 


(kj + h)(ki + kj ) 


e iki(x+y) e 2i{ki+kj)x . 


(26) 


The final result which is presented by Ecp(|2(]) is in agreement with results obtained in Ref.Q. The method used 
for deriving Eq.(^6j) enables us to investigate the sensitivity of the method for small data-errors by using Lyapunov 
exponents. 

Before we proceed, we place a few remarks. Firstly, it is interesting that large classes of solutions of the Marchenko 
equation can be obtained by solving a linear recursion relationship. The instabilities in the obtained solutions can be 
explained by the process of repeated iteration which has to be applied to compute the full kernel. This also implies 
that the stability of the obtained solutions depends strongly on the the numerical value of the parameters ki and Ai. 
The second remark is that from a conceptual point of view there is no difference in solving inverse scattering problems 
or solving nonlinear evolution equations [ p~0| . It is shown this publication that both integrable nonlinear evolution 
equations can be linearized using a similar approach as used in this paper. This explains why the inverse scattering 
method works so well in the field of integrability. 


3 Stability estimates for the Marchenko equation 


In the previous section solutions of the Marchenko equation are discussed. We have shown that we can compute 
solutions of the Marchenko equation by solving the linear recursion relation ( f25| ) . Solutions of the Marchenko equation 
can be unstable due to process of repeated iteration which is used to compute the expansion coefficients of the obtained 
solutions. The fact that the Marchenko equation can be solved by using a recursion relationship enables us to perform 
the stability analysis of the inverse scattering problem by using Lyapunov exponents. At a first sight the bookkeeping 
for the process of error propagation of the coefficients B(ni,ni,n 2 ,h 2 , ■ ■ ■ , tin, tin) is quite tedious, but the analysis 
can be simplified by considering only the coefficients which are obtained repeatedly inserting in the same recursion 
relation of Eq.@: 


B(m + 2,1,0,0, ••• , 0,0) = A\ 
B(0, 0 ,772 + 2,1, • • • ,0,0) =A 2 


fl(m, 1,0,0,-•• ,0,0) 
ifa [IV + 1] 

-8(0,0,772,1,-" ,0,0) 
ik 2 [N + 1] 


(27) 


8(0,0, 0,0, • • • , tin + 2,1) = A n 


8 (0,0,0,0, • • • ,m y,l) 
ik N [N+ 1] 
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By discussing the error-propagation in Eq. dH) instead of Eq. (f25[) has the advantage that we can investigate the 
stability of every pair (A,;, ki) in the data-spectrum. Suppose that all the parameters Ai and ki in the data lead 
to coefficients B in Eq.(^) which are stable for error-propagation. This implies that automatically all the other 
coefficients B (which are computed by solving Eq.(|25|)) are also stable for error-propagation. If there are pairs (Ai, ki) 
leading to instabilities, we can simply remove them from the data-set and hence a stable inverse scattering problem 
remains. Before computing the Lyapunov exponents, it is convenient to introduce a more efficient formulation. We 

(n) 

firstly define coefficients B\ in the following way: 


B[ n) =B(0,0,0,0, + ,0,0) (28) 

Due to imperfections in the data-set, we have to investigate the stability of the iteration series ( |25| ) under the variations 
Ai —► Ai + A Ai, ki —> ki + Aki and B^ —► B^ + A B^ respectively. If we insert these perturbations in Eq.(p5|), we 
obtain the following result for the perturbation AB \ n> : 


A B[ n) = Ci[AB\ n) } 

where the operator Ci [ • ] which expresses the error made per iteration is given by: 

, Ai B^AAi-AiB^iN + ljAki 

c>[ ' 1 = WTTE m - iNTWl - 


(29) 


(30) 


In principle, after n iterations, we obtain an error which yields: 

A Bl n+2) = C?[AAi] 

If the number of iterations n is large we can re-express Eq. (|3l| ) in the following from: 

A s| n+2) = (t\t) n av AAi 

where the averaged growth of the error per iteration (t l M ) av is given by: 


(t]Vl)av — 


Ej 1 |A^ W) | 

M 


Eq.p2) can be reformulated as: 


A B (n) = e " |og <‘»>-«A Ai = e<A4, 


(31) 


(32) 


(33) 


(34) 


where 


K = !og 


e^iab^h 


(35) 


We have obtained a closed expression for the Lyapunov exponent X l n which describes the average growth of errors per 
iteration. Although the the expression for the Lyapunov exponents may look complicated, the method can easily be 
implemented numerically. If the Lyapunov exponent is positive the errors in the coefficients grow. If the Lyapunov 
negative the errors damp out. In the latter case the inverse scattering problem is not sensitive for small errors in the 
data. In principle, we have obtained a practical approach to stabilize the inverse scattering problem since each of the 
Lyapunov exponents associated with Eq.(35) corresponds to a single Fourier component of the data-spectrum. By 
removing the Fourier components associated with a positive Lyapunov exponent, we can find a stable sub-set of the 
data-spectrum. 
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We illustrate this principle numerically for the case that that we have only one Fourier component. If we assume 
that there are only uncertainties in A B, we obtain that the Lyapunov exponent is given by: 

Al =„ 1 ™ 0 A «= 1 °g(^) ( 36 ) 

In Fig.l, the result of the reconstruction is presented for d = 1 and 0 = 10. It follows from Eq. (p6[) that for this choice 
of the coefficients d and 0, the Lyapunov exponent is negative and the inverse scattering problem is stable. If the 
experiment is repeated with the same value of 0 but with d = 0.99, the curve corresponding to the reconstruction can 
not be distinguished from the unperturbed construction. 

This situation changes if we choose 0 = 0.49 and keep d equal to unity. The reconstruction is for this case presented 
by the solid line in Fig.2. In this case the reconstruction is according to Eq.(^) sensitive to small changes in d. Indeed, 
if we solve the inverse scattering problem for d = 0.99, we find that the reconstruction given by the broken curve in 
Fig.2 differs from the unperturbed solution. 

This simple example indicates that the ratio between d and 0 determines whether the inverse scattering problem 
is stable or not. The results given in Eq.(^) indicate that a similar result holds in more general cases and forms in 
principle a tool to investigate which part of the data leads to stable inversion results. In principle if 0 —> 0, instabilities 
in the obtained solutions can be expected. However in contrast to the results presented in Ref.j^], we show in this 
paper that there is a delicate balance between the amplitude d and frequency 0 of each Fourier component to provide 
stable inversion results. 


4 Conclusions 

In this paper, we have presented a method to solve the one-dimensional Marchenko equation for finite data-sets and 
to investigate the stability of the obtained solutions. We have established a recursion relationship between the Fourier 
expansion coefficients of the data and the Fourier expansion coefficients of the kernel. We have also shown that we 
can estimate the stability of the obtained solutions from the data-uncertainties by using Lyapunov exponents. 

We want to remark that investigating the stability of the Marchenko has become an easy task from computational 
point of view since we have presented the Lyapunov exponents in an analytically closed form. A practical approach to 
find a stabilizing transformation for a potentially unstable data-set consists of simply removing all the unstable pairs 
( k n ,A n ) from the data. The low frequency components in the spectrum are responsible for the instabilities. Instead 
of bandpass filtering the data-set (to remove the low-frequency components of the data), the method presented in this 
publication enables us to identify the pairs (fc„, A n ), and their removal leads to a more accurate regularization method 
for the inverse scattering problem. 

As a concluding remark we have a renewed look on the examples of data-contaminations given in Ref. [|j . The first 
example discussed in this paper is a small DC-error on the data-set which leads unstable inversion results. It is shown 
that a DC-error corresponds with a Fourier component with vanishing k. As a result of the discussion given in this 
paper, it can be understood easily that a vanishing k leads to a positive Lyapunov exponent. The other examples 
discussed in Ref. Q deal with an amplitude error and a timing error. It is shown for both examples that if the data-set 
contains Fourier components with small k, that inverse scattering problem is unstable. The results obtained in this 
paper are in agreement with Ref. (Q , but additionally quantitative results are given. 
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Captions for Figures 

Figure 1: Kernel of K(x,x) is case of one Fourier component as given by Eq. (0) for d = 1 and f3 = 10 (solid line) 
and for d = 0.99 and (3 = 10 (broken curve). In the stable case both curves can not be distinguished. 

Figure 2: Kernel of K(x,x) is case of one Fourier component as given by Eq.(^d[) for d = 1 and f3 = .49 (solid line) 
and for d = 0.99 and j3 = 0.49 (broken curve). Both curves are scaled with a factor 10 -15 . 
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